### Replication script for Appendix A5 (robustness of correlational models) ###
library(tidyverse) # 1.3.2
library(fixest) # 0.10.4

# Load in
load("lb_clean.RData")
load("district_covars.RData")
load("district_covars_raw_ipums_codes.RData")

# Table A.18
# Combine survey data with district covars
lb_all <- left_join(lb_clean, district_covars)

presence_only <- lb_all %>%
  filter(crim_presence == 1)

m1 <- feols(crim_presence ~ confidence_government + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_all)
m2 <- feols(crim_presence ~ local_gov_responsive + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_all)
m3 <- feols(crim_presence ~ police_corrupt + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_all)
m4 <- feols(crim_gov ~ confidence_government + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = presence_only)
m5 <- feols(crim_gov ~ local_gov_responsive + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = presence_only)
m6 <- feols(crim_gov ~ police_corrupt + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = presence_only)

etable(m1, m2, m3, m4, m5, m6,
       order = c("Confident in government", "Local gov. is responsive", "Police are corrupt"),
       digits = "r3",
       cluster = ~ district_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(confidence_government = "Confident in government",
                local_gov_responsive = "Local gov. is responsive",
                police_corrupt = "Police are corrupt",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs"),
       tex = T)

# Table A.19
m1 <- feols(crim_gov ~ coercive_presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = presence_only)
m2 <- feols(crim_gov ~ presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = presence_only)
m3 <- feols(crim_gov ~ coercive_presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = lb_all)
m4 <- feols(crim_gov ~ presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = lb_all)

etable(m1, m2, m3, m4,
       order = c("State presence \\(coercive\\)", "State presence \\(all\\)"),
       digits = "r3",
       cluster = ~ district_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(coercive_presence_pca = "State presence (coercive)",
                presence_pca = "State presence (all)",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs",
                road_density = "Road density",
                log_population = "ln(Population)",
                pop_density = "Population density",
                luminosity_pc = "Luminosity per capita"),
       tex = T)

### Table A.20
m1 <- feols(crim_gov ~ confidence_government + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_all)
m2 <- feols(crim_gov ~ local_gov_responsive + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_all)
m3 <- feols(crim_gov ~ police_corrupt + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_all)

etable(m1, m2, m3, 
       order = c("Confident in government", "Local gov. is responsive", "Police are corrupt"),
       digits = "r3",
       cluster = ~ district_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(confidence_government = "Confident in government",
                local_gov_responsive = "Local gov. is responsive",
                police_corrupt = "Police are corrupt",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs"),
       tex = T)

### Table A.21
lb_subset <- filter(lb_all, !COUNTRY %in% c("Colombia", "Venezuela"))
presence_only_subset <- lb_subset %>%
  filter(crim_presence == 1)

m1 <- feols(crim_presence ~ confidence_government + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_subset)
m2 <- feols(crim_presence ~ local_gov_responsive + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_subset)
m3 <- feols(crim_presence ~ police_corrupt + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_subset)
m4 <- feols(crim_gov ~ confidence_government + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = presence_only_subset)
m5 <- feols(crim_gov ~ local_gov_responsive + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = presence_only_subset)
m6 <- feols(crim_gov ~ police_corrupt + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = presence_only_subset)


etable(m1, m2, m3, m4, m5, m6,
       order = c("Confident in government", "Local gov. is responsive", "Police are corrupt"),
       digits = "r3",
       cluster = ~ district_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(confidence_government = "Confident in government",
                local_gov_responsive = "Local gov. is responsive",
                police_corrupt = "Police are corrupt",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs"),
       tex = T)

# Table A.22
m1 <- feols(crim_gov ~ coercive_presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = presence_only_subset)
m2 <- feols(crim_gov ~ presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = presence_only_subset)
m3 <- feols(crim_gov ~ coercive_presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = lb_subset)
m4 <- feols(crim_gov ~ presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = lb_subset)

etable(m1, m2, m3, m4,
       order = c("State presence \\(coercive\\)", "State presence \\(all\\)"),
       digits = "r3",
       cluster = ~ district_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(coercive_presence_pca = "State presence (coercive)",
                presence_pca = "State presence (all)",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs",
                road_density = "Road density",
                log_population = "ln(Population)",
                pop_density = "Population density",
                luminosity_pc = "Luminosity per capita"),
       tex = T)

# Table A.23
m1 <- feols(crim_presence ~ confidence_government | district_code, weights = ~ wt, data = lb_all)
m2 <- feols(crim_presence ~ local_gov_responsive | district_code, weights = ~ wt, data = lb_all)
m3 <- feols(crim_presence ~ police_corrupt | district_code, weights = ~ wt, data = lb_all)
m4 <- feols(crim_gov ~ confidence_government | district_code, weights = ~ wt, data = presence_only)
m5 <- feols(crim_gov ~ local_gov_responsive | district_code, weights = ~ wt, data = presence_only)
m6 <- feols(crim_gov ~ police_corrupt | district_code, weights = ~ wt, data = presence_only)

etable(m1, m2, m3, m4, m5, m6,
       order = c("Confident in government", "Local gov. is responsive", "Police are corrupt"),
       digits = "r3",
       cluster = ~ district_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(confidence_government = "Confident in government",
                local_gov_responsive = "Local gov. is responsive",
                police_corrupt = "Police are corrupt",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs"),
       tex = T)

# Table A.24
m1 <- feols(crim_gov ~ coercive_presence_pca | idenpa, weights = ~ wt, data = presence_only)
m2 <- feols(crim_gov ~ presence_pca | idenpa, weights = ~ wt, data = presence_only)
m3 <- feols(crim_gov ~ coercive_presence_pca | idenpa, weights = ~ wt, data = lb_all)
m4 <- feols(crim_gov ~ presence_pca | idenpa, weights = ~ wt, data = lb_all)

etable(m1, m2, m3, m4,
       order = c("State presence \\(coercive\\)", "State presence \\(all\\)"),
       digits = "r3",
       cluster = ~ district_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(coercive_presence_pca = "State presence (coercive)",
                presence_pca = "State presence (all)",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs",
                road_density = "Road density",
                log_population = "ln(Population)",
                pop_density = "Population density",
                luminosity_pc = "Luminosity per capita"),
       tex = T)

# Table A.25
m1 <- feols(crim_presence ~ confidence_government + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs, weights = ~ wt, data = lb_all %>% drop_na(district_code))
m2 <- feols(crim_presence ~ local_gov_responsive + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs, weights = ~ wt, data = lb_all %>% drop_na(district_code))
m3 <- feols(crim_presence ~ police_corrupt + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs, weights = ~ wt, data = lb_all %>% drop_na(district_code))
m4 <- feols(crim_gov ~ confidence_government + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs, weights = ~ wt, data = presence_only %>% drop_na(district_code))
m5 <- feols(crim_gov ~ local_gov_responsive + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs, weights = ~ wt, data = presence_only %>% drop_na(district_code))
m6 <- feols(crim_gov ~ police_corrupt + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs, weights = ~ wt, data = presence_only %>% drop_na(district_code))

etable(m1, m2, m3, m4, m5, m6,
       order = c("Intercept", "Confident in government", "Local gov. is responsive", "Police are corrupt"),
       digits = "r3",
       cluster = ~ district_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(confidence_government = "Confident in government",
                local_gov_responsive = "Local gov. is responsive",
                police_corrupt = "Police are corrupt",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs"),
       tex = T)

# Table A.26
m1 <- feols(crim_gov ~ coercive_presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc, weights = ~ wt, data = presence_only)
m2 <- feols(crim_gov ~ presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc, weights = ~ wt, data = presence_only)
m3 <- feols(crim_gov ~ coercive_presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc, weights = ~ wt, data = lb_all)
m4 <- feols(crim_gov ~ presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc, weights = ~ wt, data = lb_all)

etable(m1, m2, m3, m4,
       order = c("Intercept", "State presence \\(coercive\\)", "State presence \\(all\\)"),
       digits = "r3",
       cluster = ~ district_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(coercive_presence_pca = "State presence (coercive)",
                presence_pca = "State presence (all)",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs",
                road_density = "Road density",
                log_population = "ln(Population)",
                pop_density = "Population density",
                luminosity_pc = "Luminosity per capita"),
       tex = T)

# Table A.27
lb_subset_2 <- filter(lb_all, !COUNTRY %in% c("El Salvador", "Guatemala", "Honduras", "Costa Rica", "Dominican Rep.", "Nicaragua", "Panama"))
presence_only_subset_2 <- lb_subset_2 %>%
  filter(crim_presence == 1)

m1 <- feols(crim_presence ~ confidence_government + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_subset_2)
m2 <- feols(crim_presence ~ local_gov_responsive + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_subset_2)
m3 <- feols(crim_presence ~ police_corrupt + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_subset_2)
m4 <- feols(crim_gov ~ confidence_government + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = presence_only_subset_2)
m5 <- feols(crim_gov ~ local_gov_responsive + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = presence_only_subset_2)
m6 <- feols(crim_gov ~ police_corrupt + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = presence_only_subset_2)


etable(m1, m2, m3, m4, m5, m6,
       order = c("Confident in government", "Local gov. is responsive", "Police are corrupt"),
       digits = "r3",
       cluster = ~ district_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(confidence_government = "Confident in government",
                local_gov_responsive = "Local gov. is responsive",
                police_corrupt = "Police are corrupt",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs"),
       tex = T)

# Table A.28
m1 <- feols(crim_gov ~ coercive_presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = presence_only_subset_2)
m2 <- feols(crim_gov ~ presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = presence_only_subset_2)
m3 <- feols(crim_gov ~ coercive_presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = lb_subset_2)
m4 <- feols(crim_gov ~ presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = lb_subset_2)

etable(m1, m2, m3, m4,
       order = c("State presence \\(coercive\\)", "State presence \\(all\\)"),
       digits = "r3",
       cluster = ~ district_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(coercive_presence_pca = "State presence (coercive)",
                presence_pca = "State presence (all)",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs",
                road_density = "Road density",
                log_population = "ln(Population)",
                pop_density = "Population density",
                luminosity_pc = "Luminosity per capita"),
       tex = T)

# Table A.29
lb_subset_3 <- filter(lb_all, !COUNTRY %in% c("Brazil"))
presence_only_subset_3 <- lb_subset_3 %>%
  filter(crim_presence == 1)

m1 <- feols(crim_presence ~ confidence_government + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_subset_3)
m2 <- feols(crim_presence ~ local_gov_responsive + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_subset_3)
m3 <- feols(crim_presence ~ police_corrupt + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_subset_3)
m4 <- feols(crim_gov ~ confidence_government + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = presence_only_subset_3)
m5 <- feols(crim_gov ~ local_gov_responsive + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = presence_only_subset_3)
m6 <- feols(crim_gov ~ police_corrupt + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = presence_only_subset_3)

etable(m1, m2, m3, m4, m5, m6,
       order = c("Confident in government", "Local gov. is responsive", "Police are corrupt"),
       digits = "r3",
       cluster = ~ district_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(confidence_government = "Confident in government",
                local_gov_responsive = "Local gov. is responsive",
                police_corrupt = "Police are corrupt",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs"),
       tex = T)

# Table A.30
m1 <- feols(crim_gov ~ coercive_presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = presence_only_subset_3)
m2 <- feols(crim_gov ~ presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = presence_only_subset_3)
m3 <- feols(crim_gov ~ coercive_presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = lb_subset_3)
m4 <- feols(crim_gov ~ presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = lb_subset_3)

etable(m1, m2, m3, m4,
       order = c("State presence \\(coercive\\)", "State presence \\(all\\)"),
       digits = "r3",
       cluster = ~ district_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(coercive_presence_pca = "State presence (coercive)",
                presence_pca = "State presence (all)",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs",
                road_density = "Road density",
                log_population = "ln(Population)",
                pop_density = "Population density",
                luminosity_pc = "Luminosity per capita"),
       tex = T)

# Table A.31
lb_subset_4 <- filter(lb_all, population < 250000)
presence_only_subset_4 <- lb_subset_4 %>%
  filter(crim_presence == 1)

m1 <- feols(crim_presence ~ confidence_government + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_subset_4)
m2 <- feols(crim_presence ~ local_gov_responsive + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_subset_4)
m3 <- feols(crim_presence ~ police_corrupt + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_subset_4)
m4 <- feols(crim_gov ~ confidence_government + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = presence_only_subset_4)
m5 <- feols(crim_gov ~ local_gov_responsive + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = presence_only_subset_4)
m6 <- feols(crim_gov ~ police_corrupt + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = presence_only_subset_4)

etable(m1, m2, m3, m4, m5, m6,
       order = c("Confident in government", "Local gov. is responsive", "Police are corrupt"),
       digits = "r3",
       cluster = ~ district_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(confidence_government = "Confident in government",
                local_gov_responsive = "Local gov. is responsive",
                police_corrupt = "Police are corrupt",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs"),
       tex = T)

# Table A.32
m1 <- feols(crim_gov ~ coercive_presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = presence_only_subset_4)
m2 <- feols(crim_gov ~ presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = presence_only_subset_4)
m3 <- feols(crim_gov ~ coercive_presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = lb_subset_4)
m4 <- feols(crim_gov ~ presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = lb_subset_4)

etable(m1, m2, m3, m4,
       order = c("State presence \\(coercive\\)", "State presence \\(all\\)"),
       digits = "r3",
       cluster = ~ district_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(coercive_presence_pca = "State presence (coercive)",
                presence_pca = "State presence (all)",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs",
                road_density = "Road density",
                log_population = "ln(Population)",
                pop_density = "Population density",
                luminosity_pc = "Luminosity per capita"),
       tex = T)

# Table A.33
lb_all_ipums <- left_join(lb_clean, district_covars_raw_ipums_codes)

presence_only_ipums <- lb_all_ipums %>%
  filter(crim_presence == 1)

m1 <- feols(crim_presence ~ confidence_government + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | ipums_code, weights = ~ wt, data = lb_all_ipums)
m2 <- feols(crim_presence ~ local_gov_responsive + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | ipums_code, weights = ~ wt, data = lb_all_ipums)
m3 <- feols(crim_presence ~ police_corrupt + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | ipums_code, weights = ~ wt, data = lb_all_ipums)
m4 <- feols(crim_gov ~ confidence_government + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | ipums_code, weights = ~ wt, data = presence_only_ipums)
m5 <- feols(crim_gov ~ local_gov_responsive + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | ipums_code, weights = ~ wt, data = presence_only_ipums)
m6 <- feols(crim_gov ~ police_corrupt + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | ipums_code, weights = ~ wt, data = presence_only_ipums)

etable(m1, m2, m3, m4, m5, m6,
       order = c("Confident in government", "Local gov. is responsive", "Police are corrupt"),
       digits = "r3",
       cluster = ~ ipums_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(confidence_government = "Confident in government",
                local_gov_responsive = "Local gov. is responsive",
                police_corrupt = "Police are corrupt",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs"),
       tex = T)

# Table A.34
m1 <- feols(crim_gov ~ coercive_presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = presence_only_ipums)
m2 <- feols(crim_gov ~ presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = presence_only_ipums)
m3 <- feols(crim_gov ~ coercive_presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = lb_all_ipums)
m4 <- feols(crim_gov ~ presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = lb_all_ipums)

etable(m1, m2, m3, m4,
       order = c("State presence \\(coercive\\)", "State presence \\(all\\)"),
       digits = "r3",
       cluster = ~ ipums_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(coercive_presence_pca = "State presence (coercive)",
                presence_pca = "State presence (all)",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs",
                road_density = "Road density",
                log_population = "ln(Population)",
                pop_density = "Population density",
                luminosity_pc = "Luminosity per capita"),
       tex = T)


